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Summary 


A computational model is presented for the prediction of solid/liquid phase-change energy trans- 
port including the influence of free convection fluid flow in the liquid phase region. The computational 
model considers the velocity components of all non-liquid phase change material control volumes to 
be zero but fully solves the coupled mass-momentum problem within the liquid region. The ther- 
mal energy model includes the entire domain and employs an enthalpy-like model and a recently 
developed method for handling the phase-change interface non-linearity. Convergence studies are 
performed and comparisons made with experimental data for two different problem specifications. 
The convergence studies indicate that grid independence has been achieved and the comparison with 
experimental data indicates excellent quantitative prediction of the melt fraction evolution. Qualita- 
tive data is also provided in the form of velocity vector diagrams and isotherm plots for selected times 
in the evolution of both problems. The computational costs incurred are quite low by comparison 
with previous efforts on solving these problems. 
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Nomenclature 


A 

b 

c 

C 

e 

Fo 

9 
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H 
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Pr 
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u, v 
Vi 

w 

2 , 1 / 

Greek 

a 

P 

A 

7 

r 

e 

A 

v 

P 


= finite difference coefficients 
= right hand side 
= specific heat 
= transient coefficient 
= specific energy 
= Fourier modulus 
= gravitational acceleration 
= enthalpy 

= half-height of domain 
= thermal conductivity 
= normal to interface 
= pressure 
= Prandtl number 
= Rayleigh number 
= Stefan number 
= temperature 

= Cartesian velocity components 
— interface velocity 
= width of domain 
= Cartesian coordinates 


= thermal diffusivity 
= isobaxic compressibility 
= change in accompanying variable 
= nondimensional domain width 
= modified diffusion coefficient 
= half-fusion temperature range 
= latent heat of fusion 
= kinematic viscosity 
= density 
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Superscripts 


c 

— continuity 

e , w, n, .s,p 

= geographical molecule location 

u, u,p,T 

= a and v velocity, pressure, temperature 

uu, up, uu 

= equation of first variable, multiplier 

up, TT 

of second variable 

(*) 

= nondimensionai 

Subscripts 


1,2, 3 

= solid, melt, or liquid 

d 

= dynamic 
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= east 
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= fusion 


= discrete location 
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= interfacial or dummy variable 

l 

= liquid 
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= north 
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= reference 
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= solid, or south 

sp 

= specified 
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= working variable, or west 

x,y 

= x or y direction 
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Chapter 1 


Introduction 


The Stefan problem, describing energy transport within a single component, phase change ma- 
terial, is intrinsically highly non-linear. This non-linearity is due to the compatibility constraints 
imposed at the solidification/meit front and involves the energy fluxes and front propagation veloc- 
ity at the interface location. However, neither the fluxes, the interface propagation velocity, nor the 
interface location itself are known a priori. As a result of this non-linearity, analytical solutions to 
the Stefan problem are difficult, at best, and available for only a few relatively simple configura- 
tions [1,2.3, 4, 5]. For more realistic problem specifications, discrete methods are required to effect 
the solution. Traditionally, finite difference or finite element methods have been used; the suitability 
and/or methodology of the spectral method has not yet been fully established, the enthalpy method 
of Shamsundar and Sparrow is widely used [6]. 

In the application of the enthalpy model to phase change energy transport problems, the solution 
of the algebraic system of equations consumes a considerable amount of computing time. This is at- 
tributable to the fact that complete transient histories are generally required and, more significantly, 
to the fact that Gauss-Siedel iteration or an equivalent procedure has been used to solve the equa- 
tion system. The requirement for an iterative solution procedure emerges from the highly non-linear 
character of the interface compatibility constraint [7], In the enthalpy method, the interface location 
is not tracked explicitly and its precise location can be resolved only to within one mesh spacing. 
The interface non-linearity, then, appears in the enthalpy model in the form of a highly non-linear 
equation of state relating enthalpy to temperature. In addition, this non-linear behavior is highly 
concentrated within the immediate vicinity of the phase change interface, and is extremely difficult 
to accommodate within the context of a simultaneous variable procedure for the equation system [8]. 
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This difficulty is further augmented in problems for which more than one phase change boundary 
must be accommodated within the domain. 

Williams and Curry [9] presented a more implicit procedure for the solution of the algebraic 
equation system for the case of one-dimensional phase change energy transport. In their paper, 
they attest to the difficulties involved in solving phase change problems and provide a procedure for 
implementation to multiple interface problems. They achieve this through an energy distribution 
technique in the vicinity of phase change interfaces which leads to a very complex solution procedure. 
In addition, extension of their procedure to more than one space dimension does not appear possible 
[10]. Since the majority of problems of practical importance involve two or three space dimensions, 
this is a serious disadvantage of their procedure. 

The numerical solution of solid/liquid phase change problems has been considerably simplified 
and the associated costs dramatically reduced as a result of the method proposed by Schneider and 
Raw [11] for one-dimensional problems. This procedure has also been employed in a two-dimensional 
environment by Raw and Schneider [12] with comparable cost reductions to those observed for the 
one-dimensional situation. Cost reductions are typically two orders of magnitude. All of the above- 
referenced methods and applications, however, have been restricted to conduction as the mode of 
energy transport without regard for buoyancy induced free convective motion. In practice, in terres- 
trial situations, it is frequently this free convective motion itself which is the dominant mechanism 
for the thermal energy transport. As such, it is crucial that this fluid motion predictive capability 
be available in a computational scheme in order to predict energy transport with sufficient accuracy 
so as to be of practical value. 

Unlike pure substances, multiconstituent systems do not exhibit a sharp interface between solid 
and liquid phases. In fact, due to impurities (intentional or otherwise), discrete phase change rarely 
occurs in practice. The phase change behavior of such systems depends on many factors including 
the phase change environment, composition, and thermodynamic descriptions of specific phase trans- 
formations. Moreover, solidification occurs over extended temperature ranges and solid formation 
often occurs as a permeable crystalline-like matrix which coexists with the liquid phase. 

Since they need not track phase interfaces, single region formulations are well suited for treating 
the continuous transition between solid and liquid phases, as well as the evolution of latent energy 
over a finite temperature range. Such formulations are generally developed from volume averaging 
techniques based on classical mixture theory. Detailed developments of the theory are available in 
the open Literature [13]-[18], as are applications to inert systems such as dispersed oil droplets in 
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water and fluid saturated granular materials [19]-[21]. The theories have been extended to phase 
change processes [22]-[24], although treatments have been restricted to one-dimensional, conduction 
dominated conditions. 

While continuum formulations have been shown to provide realistic predictions of transport be- 
havior for conduction phase change problems, inclusion of advective components of momentum, 
energy and species transfer does not appear to have been considered. Such an extension necessitates 
consideration of multiphase region morphology, as well as relative phase velocities. While classi- 
cal theories clearly acknowledge the significance of these factors, the desire to maintain universal 
generality prohibits description beyond that of symbolic representations. Accordingly, the primary 
objective of the present work is to develop a consistent set of continuum equations for the conserva- 
tion of mass, momentum, energy, and species in a binary, solid-liquid phase change system. Emphasis 
is placed on casting the equations into forms which are amenable to dear physical interpretation, as 
well as to solution by conventional finite-difference or finite-element methods. Although achieving 
this objective must come at the expense of a loss of generality, related assumptions and constraints 
will be clearly identified and justified on the basis of physical considerations. 
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Chapter 2 


Formulation of the Mathematical 
Model 

2.1 Binary Mixture Phase Change Problems 

In the solid /liquid, phase change of pure materials a discrete interface always separates the two 
phases. Only solid exists on one side of the interfacial surface and only liquid on the other. In a 
binary constituent phase change system the interface can take two forms. The interface may behave 
like that of a pure material as described above. Frequently however, between the region of solid and 
the region of liquid, there is a finite region with liquid finely interdispersed with solid. This region 
is often described as being ; ‘mushy”. Morphologies like this may be caused by supercooled dendritic 
growth. Dendrites can grow into each other creating a solid matrix having pockets or channels of 
entrapped liquid throughout. 

If a binary system exhibits a behavior like that of a pure material, then conservation of mass 
and conservation of momentum for the liquid region are described by the familiar two-dimensional 
equations 

+ V •(/>?) = 0, (2.1) 

= -V • (pUV) - V • f x - + pb x , and (2.2) 

at ox 

= -V • (pvV) - V • f y - + pby, (2.3) 

where r x and f y are the stress vectors defined by 

T x - T xx i + T yx j and fy = Txyi + T yy j (2.4) 
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and where b x and b y are the per unit volume body force components for the x and y directions 
respectively. The energy equation for the liquid portion is 

^il = -V-(pe l V) + V-(k i WT) (2.5) 

where e\ represents the specific internal energy for the liquid and ki is the thermal conductivity 
of the liquid. Internal energy generated by viscous dissipation and compression work have been 
neglected. Conservation of energy for the solid region does not need to include convective terms 
since the velocity of the solid is assumed zero: 

— v • (fcVT). (2.6) 

at 

In a binary system which has a region of interdispersed solid and liquid, the above equations 
are not convenient. It becomes impractical to resolve all of the details of the solid/liquid structure 
in the mushy region. The mathematical formulation for this phase change situation will be derived 
assuming solid and liquid phases can coexist in a locally distributed manner within control volumes 
of differential dimensions. 


2.2 Phase Conservation of a General Scalar Quantity 

2.2.1 Introduction and Definitions 

Mass, momentum, energy and species conservation equations are derived in this section for a single 
phase within a control volume which may contain more than this one phase. Here phase refers to 
either solid or liquid while constituent refers to either a chemical element or a compound. In this 
formulation a constituents within k phases are allowed to occupy any differential volume in space 
simultaneously. 

The control volumes used in the derivations are shown in Figs. 1, 2, and 3 with discrete interfaces 
separating solid from liquid. The precise location of the phase interfaces will generally not be known 
but the volume fractions of solid and of liquid will be used to characterize the state of dispersion at 
any location. 

Some nomenclature relating to the distribution of solid and liquid must be defined before pro- 
ceeding. The mean velocity of the constituents comprising phase k relative to a fixed reference frame 
is denoted by V^. The volume fraction of constituent a in phase k, that is, the volume occupied 
by constituent a divided by the volume of phase k which contains the volume of constituent a, is 
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denoted by g % . Similarly the volume fraction of phase A\ that is the volume occupied by phase k 
divided by the volume over which gk is evaluated which may be partially occupied by other phases, 
is denoted by g^. The actual density of constituent a in phase k is denoted by This density is the 
same that would be determined from a sample of constituent a removed from phase k. The actual 

density of phase k is denoted by p^ and describes the mass per unit volume of a volume containing 

only phase k. It should be noted that the sum over all phases of g * gives the result 

= !• ‘ ( 2 - 7 ) 

k 

The partial density of constituent a in phase k, denoted by describes the mass of constituent a 
contained by a volume occupied only by phase k, divided by that volume. This volume may contain 
other constituents than just constituent a. A mathematical expression relates partial density to 
actual density: 

Pi = 9 ipi- (2-3) 

The partial density of phase k, denoted by pk, describes the mass of phase A: in a volume which may 
contain phases other than phase k , divided by that volume. Mathematically pk is defined by 

Pk = 9kPk • (2.9) 


In a manner analogous to Dalton’s rule of partial pressures, here the sum of all partial densities at 
a location represents the actual density 

P = (2-10) 

k 


where p denotes the actual density. 

The mass fraction of constituent a in 
a in phase k, by 


phase k , /£, is related to the partial density of constituent 



( 2 - 11 ) 


This mass fraction represents the mass of constituent a in a volume entirely occupied by phase k 
divided by the mass of the entire volume. The mass fraction of phase A:, /*, is defined as 


fk = 


Pk 

ZkP'k 


( 2 . 12 ) 


the mass of phase k in a volume partially occupied by other phases divided by the mass of the total 
volume. 

The conservation of a general quantity 4> for a single phase k , <j>k, is derived here in a Cartesian 
coordinate framework for a two-dimensional control volume containing a mixture of phases. The 
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Figure 2.1: 

control volume contains two phases with the total interfacial surface between the phases denoted by 
.4/ as shown in Figure 1. For convenience the bottom left corner of the control volume is located 
at (zq,2/o). The total volume AV~ includes both solid and liquid portions and is the product of the 
width Ax and the height Ay since unit depth is considered. 

The general conservation principle is given in words by 

Time rate of change of 
<pk in volume occupied = 
by phase k 

{T 1 ) 

+ Net flow into phase k + 
across interfacial surface 
(T3) 

These four terms, the accumulation term (TT), the control face flow term (T! 2), the interfacial flow 
term (T 3), and the source term (T 4) will be treated individually in the subsections that follow. 
Expressions will be determined for a control volume of differential dimensions and the limiting result 
per unit volume will be presented as the control volume dimensions are shrunk to zero. 


Net flow of <j> across 
phase k portion of control 
faces 

(2.13) 

Source of <j>k in 
phase k volume. 

(T4) 
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2.2.2 The Accumulation Term 


The symbolic representation of the accumulation term (Tl) is given by 

Q * 

Tl = -r- I (pWkW 

dt Jv k 


(2.14) 



r 


& 


r'oN 


Figure 2.2: 

where Y\ is the portion of the control volume occupied by phase fc. In Cartesian coordinates this 
integration could be viewed as an integration between two curves, say Xi(y) and x 2 (y), or i/i(x) and 
02 ( 2 ), as shown in Figure 2. The integral part of the accumulation term can then be written in the 
forms 


t /*yG+Ay rx 2 (y) 

/ {pk4>k)dV= / / {p k 4>k)dxdy 


'yo Jx\(y) 

rxo+tix ry 2 (x) 


rxoi-Lix r 

Jx 0 Jy 


yi(*) 


(Pk4>k)dydx 


(2.15) 


Performing a Taylor series expansion of the integrand about ( xo,yo ) gives 

fyo+Ay rx 2 (y) [* Q 

[Pk<Pk)dV = I I |(Pfc*fc)|„ + 

JY k Jyo 

d 


r yyo+Ay rx 2 (y) r Q 

/ ( pk4>k)dY= / / (Pk<Pk)\ 0 + -z~(pk<t>k) \ 0 (x - x 0 ) 

JVk •Jyo Jx\(y) L 


+ fy Mk) \o ^ ” yo ^ + *** 


dxdy. 


(2.16) 


By bringing constants out of the integral and adjusting the limits, the above expression may be 
simplified to 

Q /-yo+Ay fX 2 (y) , v d 

( Pk4>k ) In / 

'yo ^(v) 


r d ryo+Ay /-X2(y) Q 

/ (/i fc 0 fc )rfy= (/?jt0fc)lo Ay fc + T“(^^)!o / / (2 ~ lo 

vX 7yo JrWy) 
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rxa+Ax ry 2(2) 

/ / {y - yo)dydx + ... 

J zn Jv\(x\ 


(2.17) 


/•xo+^x ry 2 (x) 

Jyi(x) 

It can be shown that ail terms except the first one become zero when Eq. 2.14 is divided by AY = 
Ax Ay and the limit is taken as Ax, Ay 0. The second term is considered below, although similar 
treatment applies to the third and all higher order terms. 

A new variable representing the local range of x-direction integration is introduced by the defi- 
nition 


Sx(y) = x 2 (y) - xi(y) 
Incorporating this definition we have that 

3 fyo+Ay fX2(y) 


(2.18) 


Q ry o+^y rx2{y) 

-r-(pk<P k) lo / / ( x ~ x o )dxdy 

UX Jy Q Jx\ (y) 


Q ry o+Ay 

= ^(«*»)ioy 


*i(y) 

[ ^6x 2 (y) + (n(y) - x Q )Sx(y) 


dy. 


(2.19) 


The terms Sx(y) and (xi (y) - xo) are both of the same order of size as Ax and can be expressed as 
a fraction of the control volume width Ax in the form 


Sx(y) = Ri(y)Ax and x x (y) - xq = # 2 ( 2 /) Ax. 
With substitution of these definitions the above equation becomes 

Q pyo+Ay ( , x2(y) 


(2.20) 


dx 


ryoT&y fX 2 \y) 

{pk<Pk) lo / / ( x - *0 )dxdy 

Jyo Jxi(y) 


= 5j(**t)lo l 

It is convenient to define 


yo 


rjRliy) Ax 2 + R l (y)R 2 (y)Ax 2 


dy. 


1 


f(y) = ^Rl(y) + Ri(y)R2(y) 


( 2 . 21 ) 


( 2 . 22 ) 


(2.23) 


such that the integrand is now Ax 2 /(y). The mean value theorem for integrals states that if f{y) is 
continuous in [y 0 , I/O + -Ay], then there exists some point £ in (yo, yo + Ay) such that 

l>ya + Ay 

'yo 

Using this result we can now write 

A(pfc*fc) | 0 Ax 2 Ay | y ° +Ay f(y)dy = j^(Pk<t>k) | 0 Ax 2 A »/(0; Vo < f < Vo + Ay. (2.24) 


pyo + Ay 

/ f(y)dy = A y/(0- 

Jva 
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The range of R x and R 2 is 0 < R x < 1 and 0 < Rq < 2. Further, both R x and R 2 , and hence /(£), 
are independent of Ax and A y. If this term is now divided by A V~ = AxA y and the limit is taken 
as Ax, Ay — 0, the result is 


hm — ( p k <Pk ) 
ax— o ox 
Ay— o 


Ax-Ay —(p k 4> k ) Lf(0 lim Ax = 0. 

0 Ax Ay ^ 


(2.25) 


In a similar manner all terms in the Taylor series except the first reduce to zero when this limit 
is taken. Thus the only contribution from the accumulation term, Tl, is 


Tl = £ t (pM l„|£. (2.26) 

The volume ratio of the liquid volume to the total volume as defined earlier is yj. = so that the 
accumulation term is finally given on a per unit volume basis by 


Tl - -Q t {gkPk<t>k) | 0 - 


(2.27) 


2.2.3 The Control Face Flow Term 



Figure 2.3: 

The term representing the flow across the portion of the control face occupied by phase k, T 2, is 
shown by Figure 3 to have four contributing components, one at each of the control faces. In symbolic 
form the net flow of 0 is given by 

T2 = J x — Jx+&X + jy ~ Jy+t\y (2.28) 
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Considering the 2 -component, a Taylor series expansion can be performed about xq giving 

d,J 2 : t 1 * 2 ^ Jx | 

J x — J x +Ax = Jx “ Jx + I 2 > X '*0 * 


(2.29) 


Note that the flows considered here are only those associated with phase k. The effective portion of 
the control faces for which there is a flow of <pk can be given by the ratio R times the actual area of 
the respective control faces. For the x- component flux the effective area is then R x • 1 ■ A y so that 
the flow itself becomes 


J x — l AyR x j x 


(2.30) 


where j x is the average flux over the actual area of the control volume surface occupied by phase k. 
The net ^-component flux term is then given by 

Jx “ Jx+Ax — ~~ Ax Ay-^{R X ] X ) lo ~~ 2 !^* r Q x 2 lo + * ‘ * (2.31) 

Dividing the x-component flow term by At" = AxA y, and taking the limit as Ax, Ay — 0 
reduces the term to 

lim ±(J X -J X+ Ar ) = A {RxJx ). (2.32) 

Az— 0-Y- OX 

Ay—*Q 

The '{/-component flow term similarly reduces to 

to o i(4-W = |;(^») (233) 

0 

The control face flux term (T2) for phase k is then 

T2 = - A{R xjx ) | Q - -j^(RyJy) ! 0 (2-34) 

where R x and R y are the effective ratios of areas for phase k flows at the control faces, and j x and 
j y are area averaged fluxes of the actual areas occupied by phase k. 


2.2.4 The Interfacial Flow Term 

The flow of 4>k across the phase interface surface is defined by 

T3 = [ j ki ■ a dA (2.35) 

Ja, 

where j k , is a flux across the interface and n is the normal to the interface directed into phase k. 
Beyond this representation, the interfacial flux term will not be specifically evaluated because later 
in the development this term is cancelled by a similar term associated with another phase. For 
convenience this term is denoted by 

T3 = Ii . (2.36) 
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2.2.5 The Source Term 


The source term is symbolically represented by 

T4 = / S k dV (2.37) 

•M 

where Sk 16 a per unit volume source variable. The only specific form for Sk that will be considered 
is the body force in the momentum conservation equation. By a procedure similar to the one used 
to reduce the control face flow term, the source term may be reduced to 

T4 = (g k Sk) | 0 (2.38) 

2.2.6 Composed Equation for Phase Conservation of a General Scalar Quantity 

The terms have been derived for the conservation equation of a general scalar quantity (p for a single 
phase k , fa, for a control volume of differential size containing a mixture of phases. Substituting 
these terms (T1 — T4) into Eq. 2.13 gives 

’O^idkPk'Pk) !q = ~ ~Q~^{Rxjx) lo — Qy ( RyJy ) 4" R lo 4” (9kSk) |q- (2.39) 



fa: 



Figure 2.4: 

In order to simplify this equation, an assumption is made about the effective area ratios R x and Ry. 
Figure 4 shows two cases in which the volume fraction of phase k, g k , is equal to the effective area 
ratios R x and R y . The assumption will be made here that 

Rr = Ry= 9k (2.40) 
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since in some cases it is true and since it is reasonable. Equation 2.39 now simplifies to 

d d d 

-Q t {9kPk<Pk) lo = - fc(9kjx) !o - Q^iSkjy) lo + A + (9kSk) lo • ( 2 - 41 ) 

Because the control volume has been reduced to differential size and the balance used in its 
derivation must apply for any choice of (io, yo ), the above equation can now be written as 

r\ A A 

-Q t igkPk4>k) = —Q^iakix) - Q^idkjy) + Ii + 9kSk • (2-42) 

2.3 Phase Conservation Equations for Mass, Momentum, Energy 
and Species 

Equation 2.42 represents the conservation of a general scalar quantity <fi for a single phase fc, <t> k , 
for a control volume containing a mixture of phases. The scalar quantity <t>ky and the flux terms j x 
and j y are specified in this section for the conservation equations for mass, momentum, energy and 
species. 


2.3.1 Mass 


The conservation of mass equation can be obtained from 2.42 with <f> = 1, j x = pkUk, jy = Pk^k and 
U = M k giving 

a a Pi 

(2.43) 


■^(9k Pk) = ~j- x (9kP k u k ) - J^( 9kP k v k ) + M k 


where M k represents the contributions of mass flow into the phase k portion of the control volume 
across the phase interface A[. 


2.3.2 Momentum 

The statement of z-momentum conservation can be obtained from Eq. 2.42 with <p k = u k , j x = 
p k u 2 k ~ Tfcn + P, jy = P k u k v k - T kyx , h = G kx and S k = p k g x . That is 

d d d 

-^(g k p k u k ) = --^{g k {p k u\ - r kxx + p k )} - -j^[g k (p k u k v k - r fcvi )] + G kx + g k p k g x (2.44) 

where G k x represents the production of momentum due to movement of the phase interface, g x is 
the x-component of the gravity vector, and p k is the pressure associated with phase k. 

The classical Newtonian constitutive equations need to be modified here due to spatial variations 
in the phase volume fraction. Bennon and Incropera [25] define the dilation rate of phase k as 
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v • (g k Vk) on a per unit volume of mixture basis. This was done to maintain consistency with the 
conservation of mass equation, Eq. 2.43. In their definition of dilation rate these authors use gk^k for 
the velocity terms as they arise in the mass conservation equation. Following Bennon and Incropera, 
any occurence of a velocity component will be replaced by times the velocity component. With 
this requirement the constitutive equations for phase k are then given by 


Qk^kxx — 


1 0 
-3 v ■ (s*V*) + ^($fcUfc) 


and 


r djgkVk) djjkUk) , 
gkTkyx - Pk[ Q x T Qy ]• 


(2.45) 


The statement of (/-momentum conservation follows directly from Eqs. 2.44 and 2.45 with the x 
and y subscripts, and the velocity components permuted. The y-momentum conservation equations 
are given by 


jfi(9kpk0k) = —jfelSkiPkVkVk ~ T kx y )\ ~ ~ T k yy + Pk)} + G ky + 9kPk9 y , 


(2.46) 


9k^~kxy — Qk^kyx ~ f^k 
9k T kyy = 2 fl k 


r5(</fcVjt) d(g k Uk) 


dx 


dy 


and 


|v-(^ fc V t ) + ^(s*wjk) 


(2.4' 


2.3.3 Energy 

The conservation requirement for energy can be obtained from Eq. 2.42 with <p k — s k , j x = pkU k e k — 

jy = Pk v k^k and /, = Efc yielding 

5 , , d \ ( , 8T\ 

Jt {9kPke k ) = ~ Tx [9k [pk u k e k ~ k„^j 

Here e k is the specific internal energy of phase k, that is the internal energy of a volume containing 
only phase k divided by the mass of that volume. Also k k is the thermal conductivity associated 
with phase k and E k represents the energy flow into phase k across the phase boundary, which is 
partially due to movement of the phase interface. 


d_ 

dy 


( , dT 

9k [Pk v k e k kk Qy J J 


+ E k 


(2.48) 


2.3.4 Species 


The conservation of species equation can be obtained from Eq. 2.42 with 4> k — f k , jx — Pk u kf k ~ 
PkD a k ^-,jy = PkVkfk -PkDklfc> and J ’ = MS. Hence 

-Q^{9kpkfk) ~ ~ (skPkVkfk ~ 9kPkE) k 
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(2.49) 


- ^(akPkVkfg-gkPkDt^) 

+ Mk 

where D £ is the Fickian mass diffusion coefficient associated with the diffusion of constituent a in 
phase k across the phase boundary due to the movement of the phase interface. 

2.3.5 Summary 

Mass, momentum, energy and species conservation equations for a single phase & in a mixture of 
phases, Eqs. 2.43, 2.44, 2.48 and 2.49, are conveniently summarized in Table 1. The individ- 

ual conservation equations are formed by substituting the terms in the body of the table for the 
corresponding term in Eq. 2.42 indicated by the appropriate column heading. 

Table 2.1: Summary of Conservation Equations for a Single Phase 



2.4 Continuum Conservation Equations 

Closure of the system of conservation equations for a binary mixture solid/liquid phase change 
problem requires fewer equations than the number of solid and liquid phase equations implied in 
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Section 2.3. By adding Like solid and liquid phase conservation equations together, the number of 
equations to be solved is reduced. The system of conservation equations is closed by relationships for 
phase mass fraction fk and composition f° which are discussed in Section 2.4.2, Material Properties. 
The continuum equations formed by summing over the phases are simplified in this section by 
invoking the definitions of fractional quantities, Eqs. 2.8 - 2.12, and by introducing physically 

realistic assumptions for a solid/liquid phase change system. 


2.4.1 Mass 


Summing the mass conservation equations for a phase fc, Eq. 2.43, over all phases results in 


dp 

dt 


d , . d 
■ Tx (P») - Ty^ 


(2.50) 


where p , u and v are the mass averaged quantities of density, u-velocity, u-velocity as defined by 
Eq. 2.10 and by 

V= - ^ p k V k (2.51) 

P T 

where V is the mass averaged velocity vector. The sum of the terms representing the flows of mass 
across the phase boundaries is zero (£* \'h = 0), since continuity requires that these flows, being 
internal to the control volume, occur at each other’s expense. 


Momentum 

The statement of x-momentum conservation for a mixture is obtained by summing the x-momentum 
conservation equation, Eq. 2.44, over each phase. This gives 



+ ^2 G kx + 5Z 9kPk9x (2.52) 

k k 

This equation may be simplified by using the definitions for average fluid quantities and by 
introducing new variables. By the definition of mass- averaged velocity, Eq. 2.51, the term on the 
left side of Eq. 2.52 becomes §- t (pu). The advective fluxes on the right side of the equation can be 
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decomposed into mean and relative components for each phase. The bracketed quantity of the first 
term on the right side may be decomposed to give 


Z9kPk»l = Pu 2 + YL^kPk^k ~ 


(2.53) 


by using the definition for mixture density, p, and mass-averaged velocity, V . The fourth term on 
the right side may be similarly decomposed giving 


Y^9kPkUk»k = puv + J2g kPk (u k - u)(v k - v) . 


(2.54) 


Often the relative components in Eqs. 2.53 and 2.54 may be neglected depending on the charac- 
teristics of the problem. The pressure term may be simplified by defining 

p = Yl gkPk ( 2 - 55 ) 

k 

By using the definition of mixture density p, the gravity term in Eq. 2.52 becomes 


'y ] 9kPkQx — P9x 


(2.56) 


By incorporating Eqs. 2.53, 2.54, 2.55 and 2.56, and the constitutive equation for stress Eqs. 2.45, 
Eq. 2.52 becomes 


- | ( , a v|[£( 2 «(-±v.( 5 t u ) + 

d . , d \r~*. ( (dg k v k dg k u k \ > 

- Ty ipuv) - te ? r (— + “aTJ; 


dg k u k 

dx 


d & 

- -rj^ 9kPk(uk - u) 2 ~ -q^ 3kpk{n k - u)(v k - v) 

» k J *• k * 


+ ^2 ^kx + P9x- 


(2.57) 


At this pointy physically realistic assumptions for a solid/liquid phase change system will be made 
so that Eq. 2.57 may be simplified. The convective terms Eqs. 2.53 and 2.54 will be addressed first. 
Two new variables representing the u- and v- velocity components of the liquid relative to the solid, 
are introduced here by the definitions 


u r = U{ — u 3 and v T ~ v\ — v s 


(2.58) 


By algebraic manipulation and by use of the above definitions Eqs. 2.53 and 2.54 simplify to 


- [9lPl(ui - u) 2 + g s ps{us - u 


2 + g,p 3 {u a -u) 2 = - -X- ( pfsflU ?) 



d d 

-Q^ {.9lPl( U l - - y ) + 9sps{u, - u)(v s - u)] = - — (pf s flU T V r ) 


(2.60) 


The stress terms may also be simplified for a solid/liquid system. It is assumed that the solid 
does not contribute to x-momentum through internal stresses. This requires 


Y{g s u,) = Y(g s v s ) = 0 

With this assumption the stress terms in Eq. 2.57 become 


d_ 

dx 


E (2u k (- ± V •(<?, V k + d9kUk 


d_ 

dx 


£ 

L k 


a k 


3 

fdgkVk 
V dx 


dx 


d_ 

dx 


+ 


dgk Uk 

Oy 


_ £ ( , _ d_ 

- dx (pints) - dy 


, dgivi 

~JT + 


It is further assumed that the local density gradients are negligible, that is 

V (—1=0 


V ■ V = V ■ 


E 

V jfc 


9k PkVk j _ / Pk 

p J "r vp 

Using the assumption stated in Eq. 2.61, this expression becomes 

Pi 


V • ( g k V k ) 


V-U = ^ V-{g,V) . 
P 


Similarly the partial derivatives become 


dgjui _ fn du 
dx p dx 

dgm _ p± dv_ 
dy p dy 
dgm _ fn du 
dy ~ p dy 


Substituting Eqs. 2.67 - 2.69, the stress terms, Eqs. 2.62 and 2.63 become 

2pi(-\jV-V + \ 

L \ 3 p\ f 


d_ 

dx 


2«(-3 V-(5,vi) + ^)j 


d_ 

dx 


d 


( dgivi 

dgm V 

d_ 

P 

( dv 

, du V 

dy 

Pi 

K. dx 

+ aJJ 

dy 

Pi ~ 
Pi 

Vd~x 

dy) . 


By assuming constant viscosity, the above equations sum to 


d_ 

Oy 




d 

P 

( dv 

du\ 

p 

. dy 

& T 
L Pi 

yds 

+ W. 

= Pi — 

Pi 


(2.61) 


dgiut V 
dx J . 

(2.62) 

dgiu,Y 

dy J J 

(2.63) 


(2.64) 

so that 

(2.65) 


(2.66) 


(2.67) 


(2.68) 


(2.69) 

duY 

Ox). 

(2.70) 


(2.71) 


4 d 2 u 1 d 2 v 
3 dx 2 + 3 dxdy 


+ 


d 2 u 


dy 2 

(2.72) 
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Substituting Eqs. 2.59, 2.60, 2.66 and 2.72 in Eq. 2.57 gives 
dpu dpu 2 dpuv p [4 d 2 u 


p 4 d 2 u 1 d 2 v d 2 u 

pi 3 dx 2 3 dxdy dy 2 


~ ff - h - Ty 

+ P9x + E G kx . 


(2.73) 


The fifth and sixth term on the right side of Eq. 2.73 are non-zero only in multiphase regions, where 
relative phase velocities are very small. Consequently these terms often may be neglected. 

The last term on the right side of Eq. 2.73, J2 k G kx , represents the production of x-momentum 
due to movement of the interface and shear stress at the phase interface. The former effect may be 
neglected for most physically realistic problems. Bennon and Incropera [1] proposed that Darcy’s 
law could be used to approximate the latter elfect. Darcy’s law models the force on a liquid flowing 


through a solid porous media such that 


F X — {Ql U T ) 


(2.74) 


where K x denotes the permeability in the x-direction. Knowledge of the.morphology of the solid/liquid 

interface is necessary to determine the suitability of this model as well as the coefficient I( x 

The equation for y-momentum conservation can be obtained from equation 2.73 by permuting 

the subscripts and exchanging v for u. This equation then becomes 

dpv d . . d 2 p [4 d 2 v 1 d 2 u d 2 v 

“ST = - TZ { -P vu ) - *:XP U ) + W - T 73 + 7 7777 + 777 


pi [3 dy 2 3 dxdy dx 2 


~ % ~ ~k ^ fsftVrUr) dy (5/s/rU ' ) 
+ P9y + E • 


(2.75) 


Energy 

The statement of energy conservation for a mixture is obtained by summing the energy conservation 
equation for phase k, Eq. 2.48, over all phases, giving 

J t E (9kPke k = - E E (9kpkU k e k ) ~ h 
L k J L 

- Yy E (. 9kp k v k e k ) - h ^ (2-76) 

where the mixture conductivity is defined by 

k = £ 9kk k • (2.77) 

k 
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The sum of the terms representing energy flow into a phase due to interface movement, 

zero since the exchange of energy represented by this sum is internal to the control volume. The 

advective contribution can be decomposed in a manner similar to that for the x-momentum equation 


giving 


E 9kpk( u k ~ u)(h k - h) 

(2.78) 

- k 


E 9kPk(vk - v)(h k — h) , 

k 

(2.79) 

* • 

(2.80) 


where the mixture enthalpy is defined by 

h = - ^2 9kPkhk . (2.80) 

^ k 

In the above equations, specific enthalpy h is assumed to be a suitable approximation of specific 
internal energy e. Substituting Eqs. 2.78, 2.79 and 2.80 into Eq. 2.76 yields 


I ('*) = - £ - Ty (M) + £ (* S 


k 0T 


dy ' ' dx V dxj ' dx V dy 

E gkpkiuk - u)(h k - h) 


E 9kPk(vk ~ v){h k - h) . 

. k 


(2.81) 


The temperature T can be eliminated from this expression by using the relation between temperature 
and enthalpy 

h k = / c k dT + h° k , (2.82) 

Jo 

and the identity 

(2.83) 


VT = — Vhfc = — V/i + - V(h,t - h) . 

Cfc Ck Ck 

Here c k represents the specific heat of phase k. Equation 2.81 then becomes 

d d d < JL ( IL ± — 

dt ^ dx ^ U dy ^ V dx \c k dx) dy Vc/t dy 


d \ k d 


d \k d 


- ^ E 9kPk(uk - u)(h k - h) 

L k 

- 4~ E 9kPk{ v k ~ v) (h k - h) 

°y Lt J 
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Equation 2.84 may be rewritten specifically for a solid/liquid system. The last two terms can 
be reduced by using the definitions for mass-averaged velocity, V, mixture density, p, and mixture 
enthalpy, h. This equation becomes 


d_ 

dt 


{Ph) = 


The last two terms of 


d 


d 


- 73 (puh) - 75 “ ( P vh ) + 


d ( k dh 


dx 

d_ 

^ dx 


8y 


dx \c, dx 


d 


— -7T + -7T 


dy 


d 

+ IT 


k 9 (h M 


dy Lc s dy 


- \pf,u a h s + pfiuiht - puh ] 

- [pfsVshs + pfmhi - pvh] . 

Eq. 2.85 are non-zero only in the multi-phase region. 


/ k_ dh\ 
\c a dy) 


(2.85) 


Species 


The species conservation equation for a mixture is obtained from Eq. 2.48 by summing over each 
phase. Hence 


dt 


Y (dkPk/k) 


d_ 

dx 

d_ 

dy 


Y \9kPkUkfk 


Y \9kPkVkfk 

. k 


«*» * d -£ 


( 2 . 86 ) 


The term in Eq. 2.48 representing transport of species a into phase k across the phase boundary, 
A/*., does not appear here because this transport, being internal to the control volume, is zero due 
to mass continuity. The advective contribution may be decomposed such that 


Y 9kPk u kfk = P u f° + Y 9kPk{u k - U )(/£ - f°) 

k k 

Y 9kPkVkfk = pvf a + Y 9kPk(v k ~ v)(fk - /“) 

k k 

Where the mixture concentration of species a is defined by 

r = - p Y 9kPk/k ■ 

Equation 2.86 is reduced by substitution of Eqs. 2.87 2.88 and 2.89 giving 


I - - fan - fan - £ 


d 


Y SkPkDt Y x (# ) 


L k 


d_ 

dy 

d_ 

dy 


ST~' n a df ° 

E nfk D k -J- 


d_ 

dx 


Y 9kPk(u k - «)(./“ - n 

L k 


Y 9k Pk («* “ v){fk ~ /") 


L k 


(2.87) 

( 2 . 88 ) 


(2.89) 


(2.90) 
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When Eq. 2.90 is applied to a solid/liquid system the assumption can be made that diffusion 
in the solid is negligible (D° = 0). The gradient of /£ may be decomposed into mean and relative 
components leading to 


I <^> - l + 1 w- - ■ “ d 


T v on - | tn + h ^ - n ■ 


Equation 2.90 then reduced to 


|(.n = - i o»n - 1 d-n 


d_ 

dx 

8 

dx 

8_ 

dx 

d_ 

dy 


} D ° l (/ “>] - 1 1 

} D ° 1 - r) ] - fy L' 


pd° | in 


p Da Ty ur-n 


£ - “)(/* - n 

. k 

£ - »)(/r - n 


L fe 


(2.91) 


(2.92) 


where the mixture mass diffusion coefficient is defined as 


D = 


(2.93) 


since jDf = 0. The last two terms of Eq. 2.92 may be simplified by using the definition of mass- 
averaged velocity, V, mixture density, p, and the mass fraction of a, /°\ such that 


dt 


(pn = - l (pun - 1 (pun 


d_ 

dx 

d_ 

dx 

d 


\p° | <n 


d_ 

dy l 


(p° | <n 


{ P D a -^ur-n 


d_ 
dy L 


^ or - n 


^ bwr - />/,u,/r - pun 

^ b/s^/s 0 ' - pfmft - P v f a 1 • 


(2.94) 


The last two terms in Eq. 2.94 are non-zero only in the multi-phase region. 

For a mixture containing many phases, Eq. 2.94 must be solved for all but one of the constituents, 
since one of the concentrations f a may be solved for using species conservation (i.e. Yla f a = U 
For a binary mixture then, Eq. 2.94 only needs to be solved for one constituent. 
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2.4.2 Material Properties 


Closure of the system of mixture conservation equations requires relationships for phase mass fraction 
//., phase concentration /£, phase enthalpy h k and temperature T. These relationships are obtained 
from material behavior described by the phase diagram, and from thermodynamic relationships for 
enthalpy. The required expressions are derived assuming constant specific heats, c s and C\, and 
constant latent heat, hj. Further, a linearized phase diagram is assumed. 

A eutectic binary system is first considered. Application of the proposed material relationships 
for a simpler isomorphous (non-eutectic) binary system follows. Figure 5 shows the linearized binary 
eutectic phase diagram. Several variables characterize the diagram. The slopes of the liquidus and 
solidus lines are denoted by m/ and m s respectively. The equilibrium partition ration k v is defined 
as the ratio of these slopes: 

k v = — (2.95) 

m s 

The eutectic temperature is denoted by T e and the fusion temperature as /“ — 0 is denoted by T m 
From the definition of mixture concentration (/“ = f 3 ff + lift) an d definitions of the 
variables that characterize the linearized phase diagram, expressions for the phase mass fraction 
fk and phase concentration f° may be derived. The mass fraction of solid f s for T > T e may be 
calculated from 

(2.96) 

At T = T e the solid mass fraction f s is described by 


h = 


1 - k v 


T - T] 


liq 


T-T„ 


fs= 1 


h-h 


sol 


hf 


(2.97) 


The expressions for phase concentrations are given by 

f a (2.98) 

f a (2.99) 

Calculating phase mass fraction and phase concentrations using Eqs. 2.96, 2.98 and 2.99 requires 
knowing the temperature T of the mixture. Calculations for phase enthalpy also require knowing 
temperature. Relationships between mixture enthalpy and temperature meet this requirement. The 
enthalpy for the solid phase is described by one expression for Regions i and ii in Fig. 5: 


J s 


f? = 



h s = c s T 


(2.100) 
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where it is presumed that h s = 0 and T = 0. This equation holds on the solidus as well so that 


^sol — T so \ 


( 2 . 101 ) 


where T so i is the temperature on the solidus. 

The enthalpies of the liquid phase and the mixture must be calculated according to the region 
in Fig. 5 in which the mixture belongs. The liquid phase and the mixture enthalpies for Region 
ii are discussed here. The temperature remains constant at the eutectic temperature, T e , until the 
enthalpy reaches the value 


he — ^sol 4" 


1 - 


Tsoi - Tuc 


1 - k p V T so l - Trr, 


V 


(2.102) 


Between this temperature and the liquidus, the enthalpy temperature expression involves solving 
for the roots of an expression derived in Appendix A. This expression is based on the definition of 
mixture enthalpy between the solidus and the liquidus namely 


h — f s /isol 4“ fl 


(2.103) 


Above the liquidus the mixture enthalpy becomes that of the liquid phase: 


hi = C{T + [(cj - ci)T e -f hf] 


(2.104) 


On the liquidus the above expression becomes 

*u q = ci 7] iq + f(c s - ci)T e + h f ] (2.105) 

Table 2 summarizes the enthalpy/temperature relationships for Region ii of Fig. 5. In the table, 
h 0 is defined by 

/i 0 = (c s - c,)T c = hf (2.106) 

For an isomorphous (non-eutectic) binary system or for Region i of Fig. 5 for eutectic system, the 
enthalpy /temperature expression must be modified. For low concentrations the non-eutectic system 
has a phase diagram, Fig. 6, similar to Region i of Fig. 5 for a eutectic system. The enthalpy in 
the solid phase is described by Eq. 2.100 for any system. The liquid enthalpy for the two cases 
addressed here is obtained from Eq. 2.104 by replacing T e by giving 


h { = c { T + [{c s - Cl )T so \ hf] (2.107) 

For these cases, Table 2 reduces since there is no such enthalpy range, h so i < h < A/, since h{ may 
be effectively replaced by h so [. This substitution must also be done for h 0 in Table 2 so that 

h 0 = (c 3 - c^Tsoi + hj (2.108) 
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Figure 2.5: 


Table 2.2: Enthalpy/Temperature Relationships for 
Region ii of a Binary Eutectic System 


Mixture Enthalpy h 

Temperature T 

h ^ ^sol 

A 

C| 

/&soi h ^ h e 

T e 

h c < h < h\i q 

root of AT 2 + BT + C = 0 


where A - c; — 


C = (k m - h 0 )T m + ^ 

h > /lljq 

o 

1 
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Figure 2.6: 

2.4.3 Interfacial Compatability Constraint 

The development at this point has not required that the solid/liquid interface be resolved. In some 
problems it may be desirable to explicitly resolve the interface. When the solid/liquid interface is 
resolved, the temperature and concentration distributions in the phase change material must satisfy 
interfacial mass, energy and concentration combatability constraints along the solid/liquid interface. 

The equilibrium temperature and concentration of the phase change material is determined ac- 
cording to the mixture phase diagram. The concentration of the liquid phase is given by the in- 
tersection of the temperature and liquidus lines. The concentration of the solid is specified by the 
intersection of the same temperature line and the solidus line. 

Energy and species balances applied at the interface lead to further compatability constraints. The 
first of these constraints is developed by applying an energy balance on a control volume element as 
illustrated in Fig. 7. The heat flow into the interface over time dt from the liquid phase consists of 
an advective and a convective component: 

C\rp 

q\ — - k\ dA — dt + V\e\p\dAdt (2.109) 

an i 

where n is the normal to the interface in the direction of the solid and Vj is a liquid velocity normal to 
the interface which may partly result from density differences due to the change in phase occurring 
at the interface. Similarly the heat flow leaving the interface and entering the solid is given by: 

dT 

q 3 = - k s dA — dt + V 3 e s p 3 dAdt 
on a 
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( 2 . 110 ) 




Figure 2.7: 

The difference in energy flows entering and leaving the interface results in a change of energy at the 
interface given by 

dE = pie\dAdn - p s e s dAdn (2.111) 

where dn represents the displacement of the interface occurring over time dt . Applying an energy 
balance at the interface results in 

qi-q 3 = dE (2.112) 


Substituting Eqs. 2.109- 2.111 in Eq. 2.112 gives 


, dT\ T . . dT 

+ Vieipi + k s — 

U U 1 / UTl 


V 3 t s p s 


= {piet ~ p,e,) 


dn 

dt 


(2.113) 


The development of the species interfacial compatability constraints is analogous to the above and 
results in 


- Pi A 


df£ 

dn i 


+ VJfpt + p,D, 

= (pi/r - Ps/n Y t 


VJ?p, 


(2.114) 


Conservation of mass leads to the relation 


dti 

Pi Vi - p,V a = ( pi -p s ) — . 


(2.115) 
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Substituting this equation into Eq. 2.113 gives 


'dT dT 

- k, — + k, = V 3 p 3 (-e,s) + p.e u Vi (2.116) 

on i on 3 

where e/ 3 = e; - e 3 and VJ = the speed of the interface. Similarly, substituting the mass equation. 
Eq. 2.115, into Eq. 2.114 results in 

- p,D, + P,D 3 = V 3 p 3 (-/,“) + p./gv; (2.117) 

where /£ = ff - f°. 

In Eq. 2.116, c/ 5 may be approximated by the latent heat of fusion, hj. The interfacial com- 
patability constraints for mass, momentum and energy are then 



PiVi - PsV s = (pi - 

, dn 

ps) ~di 

( 2.115) 

dT 

- *' d^ 

, dT 

+ k 3 — — — 

i dn s 

^3p3^f 4" 

(2.118) 

n 

piDi *r, 

df Q 

+ ^> -t. = 

- V sPs ft + p'fWi 

( 2.117) 


Since V s is often very small or zero, the first term on the right side of Eq. 2.118 and the first term 
on the right side of Eq. 2.117 may be neglected for many problems. 

2.4.4 Boundary Conditions 

Boundary conditions in this formulation are applied as a general convective-like boundary condition 
for any general dependent variable, cb. The form is then given by the equation 

kA - = - h* + C (2.119) 

on 

where n is the outward normal to the boundary, k $ is the conductivity of </>, and h $ is the convective 
coefficient for the quantity cf>. The variable <f> could be any of the dependent variables such as T, u, v y 
or f*. It is easily shown that through appropriate choice of the constants h ^ and C, Eq. 2.119 may 
be used to represent Dirichlet and Neumann conditions as well as the more general case as specified. 
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Appendix A 


In this appendix the relationship for mixture enthalpy is derived for a binary eutectic system in 
Region ii shown by Fig. 5, for enthalpies in the range h e < h < /iu q . The definition of mixture 
enthalpy states 

h = f s h s + fihi . (Tl) 


For a mixture whose temperature and concentration lie between the liquidus and the solidus, the 
mixture enthalpy becomes 

h = fa /l S ol + fl ^Uq • (^2) 


To proceed, we must first define each of the quantities on the right side of the above equation. The 
solid mass fraction f 3 is given in terms of the phase diagram variables by Eq. 2.96, 


f* = 


1 - k 


P L 


\ T - T Uq 

T-T m 


Further, the liquid mass fraction /; can be calculated from Eq. 2.96 since fi — 1 — f s : 


// = !-/« = 


-1 

l - k v 


Tm - T llq 


la i h 

T-T m + p 


2.96) 


(.43) 


The enthalpy of the solid on the solidus is given by Eq. 2.101 


h so i — Tsoi 


( 2.101) 


The enthalpy of the liquid on the liquidus is given by Eq. 2.105 which can be rewritten as 


hjjq — C l Tj]q -f" h Q 


(A4) 


where h 0 = ( c s - c;)T e + hf . Substituting Eqs. 2.96, 2.4.4, 2.101 and A4 into Eq. A7 and 
rearranging gives 


ci 


c e c s 

1 - kp 


T 2 + 


(c< Ca)T[jq 

1 - fcp 


1 - Zip 


ci T m 


h 


T 


+ (h- h 0 )T m + ^£2 = 0 (A5) 

1 ftp 

The root of Eq. 2.4.4 which is physically realistic gives temperature T, as a function of mixture 
enthalpy h for a mixture between the liquidus and solidus in Region ii shown in Fig. 5. This 
corresponds to the entry in Table 2 for h t < h < /iu q . 
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